Exact results for the Barabasi queuing model 
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Previous works on the queuing model introduced by Barabasi to account for the heavy tailed 
distributions of the temporal patterns found in many human activities mainly concentrate on the 
extremal dynamics case and on lists of only two items. Here we obtain exact results for the general 
case with arbitrary values of the list length L and of the degree of randomness that interpolates 
between the deterministic and purely random limits. The statistically fundamental quantities are 
extracted from the solution of master equations. From this analysis, new scaling features of the 
model are uncovered. 
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I. INTRODUCTION 

Many human activities, such as mail and e-mail ex- 
changes, library loans, stock market transactions or 
even motor activities [2| , display heavy tailed inter-event 
and waiting time distributions. To account for these 
heavy tails a priority queuing model has been pro- 
posed by Barabasi |3[, that since then stimulated an ac- 
tive field of research with potential practical applications 
(e.g., seeRefs. 0,110,1). 

Within Barabasi priority queuing model (BQM), each 
item in a list of fixed length L has a priority value. At 
each time-step, the maximal priority task is executed 
with probability p, otherwise, a randomly selected one is 
accomplished. Once a task is executed, it is substituted 
by a new one (or the same) that adopts a new randomly 
selected priority value drawn from a probability density 
function (PDF) p(x). This simple model yields power- 
law tailed distributions of inter-events times, mimicking 
the empirical histograms of many human activities. 

Besides the value of queuing models for diverse practi- 
cal questions, another issue that makes BQM attractive is 
its connection with diverse other physical problems such 
as invasion per colation [8, 9] or self-organized evolution- 
ary models [id [Til [12] , as soon as the roles of priorities 
and fitness can be identified. 

However, exact results for the BQM, both for steady (fj 
and transient [H regimes, have been obtained for the sim- 
plest instance L — 2 only. Although lists of two items 
already display the power-law decay of the distribution of 
waiting times when p approaches unity, naturally, other 
features are missed in the simplest case. Moreover, spe- 
cial attention has been given to the particular, and more 
tractable case, of extremal dynamics whenp — > 1 [l|, la |9| , 
while non-null degree of randomness (1 — p ^ 0) may also 
display interesting features. Then, in the present work, 
we tackle the BQM with arbitrary values of p and L. 

The manuscript is organized as follows. In the next sec- 
tion we show exact results for the PDFs of priorities in 
lists of arbitrary length L, by recourse to a master equa- 
tion. In Sec. Ill we obtain an approximate expression 
for the waiting time distribution. Sec. IV deals with ex- 



act results for "avalanches" which provide the time that 
higher priority tasks (above a threshold) remain in the 
list, and is also related to waiting time duration. The 
last section contains final remarks. 



II. EXACT TREATMENT 

A fundamental quantity is the probability that there 
are n tasks with priority higher than a given value x, at 
time t, P ni t(x). Its time evolution is ruled by a master 
equation (ME) of the form 

Pn,t+1 — Mn,n+lPn+l,t+Mn,nPn,t+Mn,n-lPn-l,t > (1) 

for n — 0,1,..., L, with the non-null elements of the 
tridiagonal matrix M given by 

A4-i,„(i) =px+ (1 -p)xn/L, 

M ntn (x) = P (l - x) + (l-p)(x(L - n) + (1 - x)n)/L, 
M n+1 , n (x) = (l-p)(l-x)(L-n)/L, (2) 

for n = 1,...,L, and additionally M\fi(x) =1 — 1, 
Mo t o(x) — x. Here we have taken p(x) = 1, however gen- 
erality can be recovered simply by redefining the thresh- 
old through x — > R(x) = J* p(x')dx' . 

Notice that the ME ([I])-© signals a biased random 
walk with reflecting boundaries at n — and n — L, set- 
ting the basis to write a continuum limit approximation. 
However, for arbitrary L, drift and diffusion coefficients 
are state dependent and the approach of biased diffu- 
sion successfully applied [5| to determine the scaling of 
the waiting time distribution, in other queuing systems 
with constant coefficients, becomes more tricky in the 
non-deterministic case p ^ 1. 

Then, let us find the exact steady solution of the ME 
([I])-© for arbitrary length L. By recursion, one gets 



LT(a + l)(l-x)' 



:Po(x), (3) 



(L - n)!r(a + n + 1)(1 - p)x r > 

for 1 < n < L, where a = pL/(l — p), and from normal- 
ization 



^P n (x)/P (x)Y 



(4) 



n=l 



2 



The distribution P n , given by Eqs. ©-(H]), can be used 
now to evaluate diverse meaningful quantities. In partic- 
ular, the PDF of the nth largest priority value can be ex- 
tracted from the condition J p n (x')dx' — X)m=o Pm(x), 
hence 

d 

Pn(x) = — ^ P m(x) ■ (5) 

Fig. Q] shows the exact PDFs of the two largest priori- 
ties in the list, pi(x) — P^(x) a,nd p2(x) = Pq{x) + P{(x), 
for L = 5 and different values of p, compared to the re- 
sults of numerical simulations of the BQM. 
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FIG. 1: PDFs of the largest (upper panel) and the second 
largest (lower panel) priority values, for L = 5, different val- 
ues of p indicated on the figure and R(x) — x. Solid lines 
correspond to exact results and symbols to numerical simu- 
lations of the BQM performed as in previous figures. In the 
insets, the average values are displayed as a function of p. 

In the fully random case p = 0, Eqs. ©-(01) yield 
P n {x) = Q(l - x) n x L ~ n , hence Pn (x) = L(£lJ)(l - 
x) n ^ 1 x L ^ n , in accord with straightforward combinato- 
rial analysis. In the opposite limit p — > 1, p\{x) gets 
closer to a unit step function at x = while P2{%) ap- 
proaches the Dirac delta function 5(x). This is expected 
since those tasks that have entered the list more recently 
and adopted priority values uniformly distributed in [0,1] 
have more chances to be chosen again, while the older 
tasks are more and more likely to remain in the list for- 
ever as p tends to 1, then the second priority value (and 
together with it the remaining ones) collapse to zero. 

For large enough L (namely, Lj (1— p) » 1), Eqs. ©- 
© lead to pi ~ H (x — 1+p) /p, where H is the Heaviside 



unit step function, and pi ~ (1 — p)(l/x 2 — l)H(x — 1 + 
p)/p 2 . In fact, finding directly the steady state solution 
of the ME HI)-©, in the limit of large L for fixed n (hence 
neglecting terms of order n/L), or also when p — * 1, one 
obtains a geometric progression that, for x > 1 — p, can 
be summed up to obtain the simple expression 

P (x) ~ (x- l+p)/p and (6) 

w ,ti«M:, for0 <„< £ . 

p n+1 x n 

(7) 

For x < 1 — p, all P n tend to vanish in the large L limit. 
Fig. [5] illustrates the performance of this approximation 
in comparison with exact results. The assumption n << 
L fails as soon as the probability that n > 0(1) becomes 
non negligible. For each x < 1—p, the exact P n is peaked 
around n~ (x — 1 + p)L/(l — p). The approximation 
becomes exact both in the limits of L — > oo and p — * 1. 
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FIG. 2: Probabilities Po(x) (for different values of p, upper 
panel) and P n (x) (for different values of n and p = 0.6, lower 
panel), at L — 100 and R(x) — x. Solid lines correspond to 
exact results, dashed lines to the large L/(l — p) approxima- 
tion. 

The PDF of all priorities x in the list, p(x) verifies 
y^, n — i Pn(x) = Lp(x). Its time evolution is given by 

p(x, t + 1) = p(x, t) + (p(x) - pp x (x,t)-(l- p)p(x, t))/L, 

that in the long-time limit leads to the relation 

pP (x) + (l-p)P(x)=R(x), (9) 
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where P(x) — p{x')dx' . 

Let us call old tasks those items whose priority has not 
been assigned at a given step. The cumulative PDF of old 
task priorities, 0(x), can be obtained from the relation 

LP(x) = R(x) + {L— l)0(x) (10) 

and, by means of Eq. ([9]) can be expressed as 

(L + p-l)R(x)-pLP (x) 



O(x) 



(L-l)(l-p) 



(11) 



In the particular case L = 2, Eqs. ©-(HI) give Po(x) — 
(1 + p)x 2 1 (1 — p + 2px) and recalling that its derivation 
was carried out for uniform p(x) but the general case is 
recovered simply through the mapping x — > R{x), then, 
Eq. ([TTj) allows to re-obtain the result of Vazquez Q, 
namely, 0(x) = (1 +p)R(x)/[l - p + 2pR(x)]. 




FIG. 3: Cumulative PDFs of old task priorities, for L = 20, 
different values of p and R(x) — x. Symbols correspond to 
numerical simulations of the BQM performed as in previous 
figures, black lines to the exact results from Eq. (|lip . 

Fig. [3] illustrates the behavior of O(x) for different val- 
ues of p and L = 20. The distribution of the bulk of 
old values for arbitrary L is qualitatively similar to that 
obtained for L = 2 in Ref. Q. 

For L > 2, however, P(x) and 0(x) are mean-field 
quantities, while more meaningful is the distribution of 
the largest old priority Oi(x) (that, of course, for L = 2 
coincides with 0{x)). It verifies 



P (x) =R(x)0 1 {x), 



(12) 



since the probability that there are no tasks above x, 
Po(x), is the product of the probability that the freshly 
assigned (new) priority value is below x times the prob- 
ability that the highest old task priority (hence also the 
remaining ones) is below x, as soon as the new priority 
value and the old ones are independent. For the partic- 
ular case p — 0, 0\(x) — R L ~ 1 (x) while in the opposite 
limit p — > 1, it tends to a unit step function at x = 0. 



For L > 2, the distribution of the second largest old 
priority 02(2) can be extracted from the identity 



P 1 = (1-R)0 1 +R(0 2 -0 1 ), 



(13) 



which comes from considering that the probability that 
there is only one task above x, Pi{x), is prob.[new > x A 
1st old < a;]) + prob.Qnew < x A 2nd old < x A 1st old 
> x\), while prob.(2nd old < x A 1st old > x)=prob.(2nd 
old < x)- prob.(2nd old < x A 1st old < x)~02{x) — 
Oi(x), since the first and second largest old values are 
not independent. Analogously, in general one has 

P n = (1 - R)(O n - O n _i) + R(O n+1 - O n ), (14) 

for 1 < n < L — 1, taking Oq — and Ol = 1, while 
Pl = (1 — -R) (1 — Ol-i)- From where the whole family of 
stationary old task distributions can be straightforwardly 
obtained. 

Exact results for Oi(x) and 02(x) are compared to the 
outcomes of numerical simulations of the BQM in Fig. [4] 
for L = 20, different values of p and R(x) = x. Observe 
that 02(x) is bounded from below by 0\{x) and more 
generally Oi(x) < 2 (x) < ... < L ^x{x). 
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FIG. 4: Cumulative PDFs of the first (upper panel) and 
second (lower panel) largest old task priorities, for L — 20, 
different values of p and R(x) = x. Symbols correspond to 
numerical simulations of the BQM performed as in previous 
figures, black lines to the exact results from Eqs. (|12[) -(U3 [ ). 
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III. WAITING TIME DISTRIBUTION 

The family of distributions of old values {O n , 1 < n < 
L — 1)} should in principle allow to compute exactly the 
distribution of waiting times P w {t) in the steady regime. 

Pw(t) can be obtained as 



1.0 



P w (r) 



dR{x)r T {x) 



(15) 



where r T (x) is the probability that a task (let us call it 
X) with freshly acquired priority value X at a given time 
t = t (once attained the steady state) is again selected 
for the first time at t = to + r. 
For r = 1, 



r 1 (x)= P 1 (x) + {l-p)/L. 



(16) 



Hence P w (l) = l/L when p = and it tends to one in 
the opposite case p — > 1. By means of the approximate 
Eq. © forPo(x), one has P W {1) ~ p+ (1 -p) ln(l -p) + 
(l-p)/L. 

The probability that, instead of X, the first old task is 
selected at to + 1 is p(l — 0\(x)) + (1 — p)/L, while the 
probability that any other old task is selected (for L > 2) 
is (1— p)/L. Given each of these L — 1 cases, for comput- 
ing 7*2(3;), one has in principle a different probability of 
selection of X at the second step (t = to + 2) which will 
be a function of 0\ and 2 - More generally, the exact 
calculation of P w (r) , for r > 1 , will require to consider a 
branching process, with L — 1 paths at each node, such 
that for L > 2, r r (x) does not factorize. This tree gen- 
eralizes the branching process considered in analogy to 
invasion percolation for L = 2 Q. 

In the first steps (up to r ~ L), the statistics will be 
conditioned by the memory of previous selections (ag- 
ing regime). This is because recently chosen tasks have 
propensity (the higher, the closer p to 1) to be chosen 
again, dominating P w at small r. 

In particular, for r = 2 one obtains 

r 2 (x) = c(l-r 1 )+pR[(p+c)0 2 +{c{L-2)-p}0 1 ], (17) 

where c — (1 —p)/L. For the special case L = 2, Oj> = 1, 
then one recovers the expression found in Rcf. |6j for 
L = 2, namely, r 2 = (1 - n)[pJ2 + (1 - p)/2]. 

At further time-steps this is a hard trail to proceed 
and the results may not be expressible in a readily man- 
ageable form. However, notice that while for small r, 
the integral of is dominated by large values of x, due 
to the propensity of such values to be re-chosen early, 
contrarily, for large enough r, r T (and hence P w ) will 
gain the main contribution from the purely random (un- 
conditioned) selection from the bulk of relatively small x 
values (as can be seen in Fig. [5] where r T {x) is displayed). 
This is expected to apply also when p — > 1 at any r. For 
such cases, one can write 

r r (a:)=i(l-ri(a:))(l-/(i)r- 2 /(a;), (18) 



o r, 


L = 3, p = 0.6 
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FIG. 5: (Color online) Integrands r T of the distribution of 
waiting times, for values of r indicated on the figure and two 
couples of parameters (L,p). Symbols correspond to numer- 
ical simulations, solid lines to exact results and dashed red 
fines to the approximate analytical expressions. 



where f(x) is the effective probability that task X is se- 
lected at some given step t > to + 1 and can be estimated 
as f(x) — pPo(x) + (1 —p)/L, as soon as P = RO\ is the 
probability that there are no tasks with priorities higher 
than x. Fig. [5] also exhibits the comparison between exact 




FIG. 6: Distributions of waiting times, for fixed size L = 10 
and different values of p indicated on the figure (upper panel) 
and for fixed p = 0.6 and increasing values of L indicated on 
the figure (lower panel), R(x) = x. Solid lines join the analyt- 
ical results from Eqs. (|15[) an d (|18|l and symbols correspond 
to numerical simulations of the BQM. Insets: rescaled plots 
of the numerical histograms, where to = l/\n(L/(L — 1 +p)). 
Dashed lines with slopes -1 (upper inset) and -3/2 (lower in- 
set) are drawn for comparison. 
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and approximated functions r T . 

In particular, for p = 0, Eq. (JTSJ) is independent of the 
choice of f(x) and it correctly yields the pure exponential 
decay P w {t) = (1 - l/L^/L for all r [3]. In the op- 
posite case p — ► 1, and using the approximation given by 
Eq. © for Po i Eq. (fT5j) leads to the asymptotic behavior 

P w (r) ~ iexp(-T/r ), (19) 

T 

where t = l/ln(L/(L- l+p)) ~ L/(l-p). 

This expression for the characteristic time To applies 
por any p. Thus, the characteristic exponential decay 
time tq is shifted to larger r when p — ► 1 as well as when 
L increases. 

Analytical predictions are compared to numerical sim- 
ulations in Fig. One observes that the approximate 
expression derived from Eq. (|15p manages to describe 
the exponential cutoff in all cases and the scaling regime 
in the limit p — ► 1, although it fails to predict the -3/2 
power-law neatly observed in numerical simulations for 
< p < 1 as L - > oo (notice in the lower panel of 
Fig. [6] the deviation for r < L, leading to a spurious 
power-law exponent -2). This is due to the fact that the 
aging regime is overlooked by this approximation. Let 
us remark that a -3/2 exponent is also found in classical 
queuing models with fluctuating length [H, Q and the re- 
turn time distribution of a random walk is at its origin. 
In view of the difficulties to find the exact expression for 
P w (t), to explain this scaling regime, we will solve next 
a closely related problem. 

IV. AVALANCHES 

Let as also consider now the events between two suc- 
cessive times when the number n of priorities above a 
given threshold x vanishes (avalanche). Avalanche du- 
ration is relevant in the present context as as soon as 
it provides the duration of intervals in which there are 
queued tasks with priorities above a threshold to be ex- 
ecuted. From the viewpoint of random walks, this is a 
first passage problem. Following the lines in Ref. [12|, let 
us define Q n ,t{x), the probability of having n values with 
priorities higher than x, given that an avalanche started 
at t = [t time units ago). Q n ,t follows the same ME 
(HI)-® as P n .t does, except for Mo,i = 0, and the initial 
condition is Qi,o = 1 — X and Q n _a = 0, Vrt > 1. Thus, 
the probability that an avalanche, relative to threshold 
x, has duration t is 

q t (x) = xQ 1)t -i(x). (20) 

Fig. [7] illustrates the scaling that comes up for any p 
at the critical threshold x = 1 — p. Exact results were 
obtained by numerical integration of the ME for Q n and 
compared to the results of numerical simulations of the 
BQM. Notice that the scaling region increases with L and 
shifts towards larger times as p — > 1 . 
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FIG. 7: Distribution of avalanche size q(t) = qt(x — 1 — 
p) for L = 100, different values of p indicated on the figure 
and R(x) = x. Solid lines join the exact values and symbols 
correspond to numerical simulations of the BQM. In the inset, 
exact results for p = 0.6 and different values of L indicated 
on the figure are displayed. Dotted straight lines with slope 
-3/2 are drawn for comparison. 

The ME of Q n can be solved analytically through di- 
verse standard methods [l3l.[T3|. Yet, in the limit of large 
L and fixed n, the ME describes a simple biased random 
walk, with an absorbing boundary at n = and proba- 
bilities to step either to the right, to the left, or remain 
still, given by to + = (1 — p)(l — x), mT — px, m° = 
1 — m + — to - , respectively. From this viewpoint, qt(x) is 
the probability that the first return to the origin occurs 
at t when the avalanche started at t = 0, while Qi,t{x) is 
the probability of reaching n = 1 at time t — 1, without 
having visited n < 0. Thus, q just differs from Qi in 
appending the last step from 1 to 0. For any n, Q n ,t(x) 
can be found by solving first the unbounded problem and 
then resorting to the reflection principle [HI, [l6| . More- 
over, if we are concerned with the asymptotic behavior, 
we can directly take advantage of the Gaussian approx- 
imation from the central limit theorem. Therefore, one 
has 

(n-l-ct) 2 (n+l-ct) 2 

Q n>t {x)~{l-x) 7 =^= , (21) 

V zna^t 

where c, and a 1 are the mean and variance of each single 
step. This readily leads to the asymptotic behaviors 

/ 1 ~ 3/2 ' if c = 

9tW ~ { ;-i/2 exp( ^ ); otherwis e, (22) 

that is, an exponential decay dominates the long-time de- 
cay in the biased cases, meanwhile, if c = to + — mT = 
(hence x = 1 — p), a power-law arises in the large L 
limit, in agreement with the results displayed in Fig. [7] 
and with the well known results for a driftless random 
walk [l6l |. In particular, there is a correspondence with 
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the random annealed Bak-Sneppen model, where the 
same scaling is observed for any K at the critical thresh- 
old x = 1 /K [l2j ■ Let us remark that in the Bak-Sneppen 
model the transition matrix for the associated ME has 
2K non-null diagonals, and a generic univoque relation 
between p and K does not emerge. However, concerning 
avalanches, the equivalence between both models arises 
for K = 1/(1 — p). Due to the threshold being an upper 
or lower bound in each case, that relation is complemen- 
tary to K = 1/p which arises by identifying ratios of 
deterministic/random sites [l|. 

V. FINAL COMMENTS 

Summarizing, we obtained analytical results for the 
BQM with queues of arbitrary length. Exact expressions 
were shown to be in agreement with the outcomes of nu- 
merical simulations of the dynamics. Progress has still 



to be made to obtain the exact waiting time distribution 
that displays different regimes between the purely expo- 
nential one (at p = 0) and the power-law decay with unit 
exponent (at p — > 1), when L — > oo. However, an approx- 
imate expression has been found that accounts for most 
of the distribution traits. Moreover, we have shown that 
avalanches, at the critical threshold x = 1 — p, constitute 
another scale- free feature of the BQM for p > 0. Be- 
sides the main applications here illustrated, the present 
results may allow to estimate many other relevant sta- 
tistical quantities of the BQM and can be extended to 
other queuing systems. Furthermore, our exact results 
set the basis to further explore the correspondence be- 
tween BQM and other related models. 
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